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Abstract To derive the convergence field from the gravitational shear 7 of the background galaxy images, the 
classical methods require a convolution of the shear to be performed over the entire sky, usually expressed thanks 
to the Fast Fourier transform (FFT). However, it is not optimal for an imperfect geometry survey. Furthermore, 
FFT implicitly uses periodic conditions that introduce errors to the reconstruction. A method has been proposed 
that relies on computation of an intermediate field u that combines the derivatives of 7 and on convolution with a 
Green kernel. In this paper, we study the wavelet Helmholtz decomposition as a new approach to reconstructing 
the dark matter mass map. We show that a link exists between the Helmholtz decomposition and the E/B mode 
separation. We introduce a new wavelet construction, that has a property that gives us more flexibility in handling 
the border problem, and we propose a new method of reconstructing the dark matter mass map in the wavelet 
space. A set of experiments based on noise- free images illustrates that this Wavelet Helmholtz decomposition 
reconstructs the borders better than all other existing methods. 
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1. Introduction 

Weak gravitational lensing provides a method of di- 
rectly mapping the distribution of dark matter in the uni- 
verse (Bartelmann & Schneider 1999; Mcllicr 1999; Van 
Waerbeke et al. 2001; Mellier 2002; Refregier 2003; Munshi 
et al. 2008; Pires et al. 2010). This method is based on 
the weak distortions that lensing induces in the images of 
background galaxies as it travels toward us through in- 
tervening structures. The induced shear is small, typically 
around an order of magnitude less than the ellipticity seen 
in background galaxies, therefore the weak lensing effect 
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cannot be measured on a single galaxy. To measure the 
shear field, it is necessary to measure ellipticities of many 
background galaxies and construct a statistical estimate 
of their systematic alignment. This measured shear field 
is directly related to the distribution of the (dark) matter 
and thus to the geometry and the dynamics of the uni- 
verse. As a consequence, weak gravitational lensing offers 
unique possibilities for probing the statistical properties 
of dark matter and dark energy in the universe. 

The mapping of the dark matter distribution has be- 
come a central topic in weak lensing since the very first 
reconstructed 2D mass maps demonstrated that you can 
see the dark components of the universe. Recently, this 
field has seen great success in reconstructing the deepest 
and largest 2D dark matter distribution in the COSMOS 
field (Massey et al. 2007). From this distribution, the na- 
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ture of dark matter can be better understood, so better 
constraints can be set on dark energy (Schrabback et al. 

2010) . This method is now widely used but, the amplitude 
of the weak lensing signal is so weak that its detection re- 
lies on the accuracy of the techniques used to analyze the 
data. Each step in the analysis has required the develop- 
ment of advanced techniques dedicated to these applica- 
tions. 

The reconstruction of the dark matter mass map from 
shear measurements is a difficult inverse problem because 
of observational effects such as noise and the complex ge- 
ometry of the field. Classical methods based on Kaiser 
& Squires (1993) are not optimal. Indeed they require 
a convolution of the observed shear with a kernel to be 
performed over the entire field, usually expressed in the 
Fourier domain, and the FFT imposes a condition of peri- 
odicity that is not true and introduces a mixing between 
E and B modes. An alternative interesting method was 
proposed in Seitz & Schneider (1996), which consists in 
convolving the shear field with a local kernel that depends 
on the shape of the domain. It was shown to be optimal in 
Lombardi & Bcrtin (1998), and it was reformulated later 
in Seitz & Schneider (2001). 

This approach is relatively expensive in computational 
resources. Much effort has been made in past years to 
properly derive a cosmic shear two-point correlation on a 
finite interval (Kilbinger et al. 2006; Schneider & Kilbinger 
2007; Fu & Kilbinger 2010). In Fu & Kilbinger (2010), a 
new statistic is proposed, maximizing the signal-to-noise 
ratio and a figure of merit based on the Fisher matrix 
of the cosmological parameters Q m and ag. However, this 
approach does not allow reconstructing a map, and can- 
not be used for higher order statistics studies such as bis- 
pectrum analysis (Van Waerbeke et al. 2001; Pires et al. 
2009b; Munshi et al. 2011), peak counting (Marian et al. 

2011) , Minkowski Functionals (Kratochvil et al. 2011), or 
wavelet peak counting (Pires et al. 2009a). 

In this paper, we study a new approach for weak lens- 
ing mass inversion in the ideal case of noise free data based 
on a Helmholtz decomposition in the wavelet domain. We 
first show that there is an explicit relation between the 
weak lensing E/B mode decomposition and the Helmholtz 
decomposition which decomposes a vector field into two 
components, one curl-free and the second divergence-free. 
In contrast to the method presented in Seitz & Schneider 
(1996), our method does not rely on a convolution with a 
kernel, and no intermediate field u derived from the shear 



field is needed. The Helmholtz decomposition is performed 
in the wavelet space, and the Fourier transform is never 
used. This gives us more flexibility to handle borders or 
to add some constraints in our solution. We introduced a 
new wavelet construction based on border-wavelets, and 
we present two kinds of wavelet Helmholtz decomposition 
algorithms. The first one imposes a global zero B modes 
constraint and the second a local zero B mode constraint 
only on the borders. 

This paper is structured as follows. In Sect. 2, we re- 
view the basis of the weak lensing formalism. In Sect. 3, 
we introduce of the Helmholtz decomposition to the weak 
lensing E/B modes decomposition, and the discretization 
is discussed in Sect. 4. Then, in Sect. 5, the Helmholtz 
decomposition is extended to the wavelet domain and, we 
present new constructions of divergence-free and curl-free 
wavelets that satisfy suitable boundary conditions. And 
finally, in Sect. 6, we do some numerical experiments for 
weak lensing and conclude in Sect. 7. In Appendix A, 
we provide the mathematical expressions for the vector 
wavelets and in Appendix B, we present the algorithm 
to perform the E mode convergence reconstruction in a 
bounded domain using a wavelet Helmholtz decomposi- 
tion. 

2. Weak lensing formalism 

2.1. The weak lensing equations 

In weak lensing surveys, the observable shear field 7 4 
can be written in terms of the lensing gravitational poten- 
tial ip(0i) (Bartelmann & Schneider 1999): 

72 = d&ip, 

where the partial derivatives c\ are with respect to 0j. 
Notice the difference in sign with respect to (Seitz & 
Schneider 1996; Lombardi & Bcrtin 1998). 

The shear 7i(6*) with i = 1, 2 is derived from the shapes 
of galaxies at positions in the image. The 71 component 
corresponds to a deformation in the horizontal/vertical 
direction and the 72 component to a deformation in the 
diagonal direction (see Fig. 1). 

The convergence n(9) can also be expressed in terms 
of the lensing potential ijj, 

*=\{ft + %)1>, (2) 
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and is related to the mass density £(#) projected along 
the line of sight by 

E(0) 



k{9) = 



(3) 



(4) 



v ' 

^crit 

where the critical mass density £ cri t is given by 

crtt ^GDolDls' 
where G is Newton's constant, c the speed of light, 
and Dos j Dol, and are respectively the angular- 

diameter distances between the observer and the galaxies, 
the observer and the lens, and the lens and the galaxies. 

2.2. Difference between the sheared galaxies and the 
shear field 7 

The deformation of a circle C by a shear field provides 
an ellipse MC and corresponds to the local linear applica- 
tion of the matrix 



(5) 



cos 2(f) sin 2(f) 
sin 2(p cos 2(f) 

where k is the convergence, <f> the angle formed by the 
larger axis of the ellipse and the {Ox) axis, and I7I the 
amplitude of the deformation. 




Figure 1. Galaxy shapes and corresponding gravitational 
shears 7 for a circular galaxy. 



In the weak lensing approximation, we only consider 
the shear 7 = |7|(cos 2(f), sin 2(f) which concentrates the 
information on the angle of deformation and the amplitude 
of deformation and defines a 2D vector field: 7 : SI — > M 2 , 
see figure 1. Although the shear is often represented by 
sheared galaxies which is a spin-2 representation, as one 
can see in the first row of Fig. 1, the shear field 7 itself is 
a usual vector field, as seen in the second row of Fig. 1, so 
we can compute its divergence and its curl: 

div7 = <9i7i + <9 2 7 2 , curl 7 = <9i7 2 - d 2 7i, (6) 

and apply the Helmholtz decomposition to it. 



2.3. The mass inversion 

The weak lensing mass inversion problem consists in 
reconstructing the projected (normalized) mass distribu- 
tion k{9) from the measured shear field ji(9). We invert 
Eq. (1) to find the lensing potential V an d then apply 
formula (2) following the classical methods based on the 
pioneering work of Kaiser & Squires (1993). In short, this 
corresponds to (Kaiser & Squires 1993; Starck et al. 2006): 



ft = A" 1 ((dj - a|) 7 i + 29x9272) 



dl - d\ 2o>id 2 

di + di 11 + W+dl 



72- 



(7) 



This method usually is called a Fourier transform and has 
been generalized to spherical harmonics in Pichon ct al. 
(2010) to deal with full-sky maps. To recover k from both 
71 and 72 there is a degeneracy when fci = fc 2 = 0. 
Therefore, the mean value of n cannot be recovered from 
the shear maps. This is known as the mass-sheet degener- 
acy. 

The most important drawback of this method is that 
it requires an observation of the shear over the entire sky. 
As a result, if the observed shear field has a finite size or 
a complex geometry, then the method can produce arti- 
facts on the reconstructed convergence distribution near 
the boundaries of the observed field. The error caused by 
the artificial periodic boundary conditions can not be de- 
creased by taking mirror conditions which remove the dis- 
continuity. 

There are local inversions that reduce the unwanted 
boundary effects (Seitz & Schneider 1996). But regardless 
of the local formula used, the reconstructed field has more 
noise than when obtained with a global inversion, and the 
inversion is unstable in the presence of uncorrelated noise. 

The convergence k can be computed directly in direct 
space (without the Fourier transform) thanks to the kernel 
integration (Seitz & Schneider 1996): 



k(6)-k = ^ J K{e,e')-j{e-e')de', 



(8) 



e' en 



where stands for the mean value of k. The kernel K 
depends on the geometry of the domain ft. For f2 = R 2 , it 
is given by 



K{<p,6) = 



®2 — ^1 — 2$i$2 



(9) 
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In Seitz & Schneider (1996), the authors propose to com- 
bine the derivatives of (7$) 



u = 



01 71 

0172 



0272 

027i 



(10) 



and then to apply the Helmholtz decomposition u = 
Vk' e ' + V x k( b ) , in order to reconstruct the convergence 

In the following, we propose to use the Helmholtz de- 
composition directly on (7$) to perform the E/B mode 
decomposition, and this can be seen as a kind of short- 
cut in comparison with the Seitz & Schneider method. 
Furthermore, we also propose to perform the Helmholtz 
decomposition in the wavelet space using a new kind of 
border- wavelet, which allows us to reduce border effects 
during the reconstruction. 

2.4. The E and B mode decomposition 

Just as a vector field can be decomposed into a gradient 
or electric (E) component, and a curl or magnetic (B) 
component, the shear field 7,(0) can be decomposed into 
two components, which for convenience we also call (E) 
and (B). The decomposition of the shear field into each 
of these components can be easily performed by noticing 
that a pure E mode can be transformed into a pure B 
mode by a clockwise rotation of the sheared galaxies by 
45°: 71 -> -72, 72 -> 71 ( scc %• !)■ 

Because the weak lensing arises from a scalar poten- 
tial (the Newtonian potential <&), it can be shown that 
weak lensing only produces E modes. As a result, the least 
square estimator of the E modes convergence field is sim- 
ply 

k e = A- 1 ((0? - 0|) 7 i + 20i0 272 ) • (11) 

The E mode deformations 7^) which satisfy (1) can be 
equivalently defined by 



7 (B) e T E = {7 : 20i027i 



2 2 ) 7 2 =O}. (12) 



On the other hand, residual systematics arising from 
imperfect correction of the instrumental PSF or telescope 
aberrations, generally generates both E and B modes. The 
presence of B modes can thus be used to test for any 
residual systematic effects in current weak lensing surveys. 
For this purpose, the following estimator for the B mode 
convergence field can be formed (Starck et al. 2006), 



and it is consistency with zero will confirm the absence of 
systematics. 

The B mode deformations 7^) derive from potentials 

tpB as 



j[ B) = 0102^3, 

-I(02-0|)V> B , 



(B) 
72 



(14) 



and form the set 



k B = A- 1 (20x0271 - (0? - 2 2 )7 2 ) 



(13) 



7 (B) G T B = {7 : (0 2 - 2 2 ) 7 i + 2010272 - 0} . (15) 



If we consider 7 defined on the infinite plan M 2 or on 
the torus T 2 , i.e. 7 e (L 2 (M 2 )) 2 or 7 e (L 2 (T 2 )) 2 , then 
T^ and Tb define orthogonal spaces for the usual scalar 
product in L 2 . This is no longer the case for domains with 
boundaries. In this case, F^ and Tb have nonzero common 
elements. 

Solutions (11) and (13) are therefore no longer valid 
for bounded domains. The gravitational shear that derives 
from the bi-harmonic function ip such that A 2 ip = pro- 
duces modes that both remain to Te and T b (Bunn et al. 
2003). Indeed, the borders cause mode-mixing or, equiv- 
alently, a leakage of E modes into B modes. Therefore, 
even a pure gravitational shear will produce both E and 
B modes in a finite field. Consequently, it will be more dif- 
ficult to confirm the absence of residual systematics with 
the presence of border effects. More details can be found 
in Crittenden et al. (2002); Schneider et al. (2002). 

In the next section, we propose to explore the use of 
the Helmholtz decomposition to reduce boundary effects 
in E and B mode decomposition. 



3. Introduction of the Helmholtz decomposition 
into the weak lensing formalism 

In this section, we introduce the relation between the 
Helmholtz decomposition and the weak lensing E/B mode 
decomposition. The Helmholtz decomposition in two or- 
thogonal components only stands for unbounded domains. 
Then the unicity of the ke reconstruction on bounded do- 
mains can only be obtained under some constraints. 
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3.1. The weak Sensing E/B mode decomposition 

On (L 2 (M 2 )) 2 , we can express the relations (1) & (14) 
thanks to the Fourier transform: 



7< E >(fc) = 
(B) (fe) = 



7 



\(k\ k 2 ) 

k\k 2 
2 ( — ^1 + ^2 



^ (E) (fc), 

v;( B )(/c). 



(16) 



It allows us to derive the convergence ke deriving from the 
E modes (the actual convergence) and the convergence k b 
deriving from the B modes (the systematic errors) simi- 
larly to (11) and (13): 



k(k) 



k E {k) 


1 


_ k B {k) _ 





kf - k{ 2kik 2 
2k 1 k 2 -kl + kl 



71 W 

72 (k) 



(17) 



where hat symbols denote the Fourier transforms. 



with: 

P7(fc) = 

®y(k) = 



l 

W 



k 2 



—k\k 2 
hk 2 kf 



7 (fe) 



1 

w 



A r 

kf hk 2 
k x k 2 



k\ 



7 (fe). 



(20) 



(21) 



From equations (17), (20), and (21), we have A K = — Ap + 
An. Then 



K = 



k E (k) 
k B (k) 



-A v 



71 (*) 

72 (k) 



■At 



71 (fc) 

72 (fe) 



(22) 



This provides the following E and B modes conver- 
gence: 



if 7 = P7 
then k = 



K E 
k b 



= -P7 + Q7. 



(23) 
(24) 



3.2. Relation between Helmholtz and E/B modes 
decompositions on an unbounded domain 

The Helmholtz decomposition theorem says that every 
vector field 7, defined everywhere in space can be decom- 
posed into a rotational part 7 div and an irrotational part 
7curio- The Helmholtz decomposition consists, for a vector 
field 7, in writing 



7 — 7div + 7curl 



d 2 ip 



( dip 

V d 2V 



(18) 



Vp with ip and p 



where 7 div0 = V x V and 7 curl0 
differentiable functions. 

This decomposition is unique for domains without 
boundaries. In the Fourier domain, it consists in decom- 
posing the vector field 7(fc) in the longitudinal direction, 
i.e. parallel to k and in the transverse direction, i.e. per- 
pendicular to k. The longitudinal component corresponds 
to the curl-free component and the transverse component 
to the divergence- free component. 

Using the following notations: P7 = 7 div0 and Q7 = 
7curio> we obtain 



3.3. Relation between Helmholtz and E/B mode 
decompositions on a bounded domain 

On the square f2 = (0, l) 2 , the Helmholtz decompo- 
sition becomes much more complicated, see (Dautray & 
Lions 1984, for further details). Indeed, the unicity of 
the Helmholtz decomposition 7 = P7 + Q7 disappears. 
However, it can be regained by considering three compo- 
nents in the orthogonal decomposition: 



7 - 7div0 +7curl0 +7ha 



(25) 
• n = on 



7 = P7 + ' 



(19) 



with div7 div0 = on the square Q and 7 d ; v0 
the borders of the square curl7 curl0 = on and 
7curio ' t = on dft, and div7 harm = curl7 harm = 
on n, where n and t are the normal and tangential 
unit vectors on dil. This new term 7h arm derives from 
a harmonic function on 0: 7 harm = Vpo = V x w 
with Apo = Aluq = 0. Alternatively we can write 
7harm = ^Po + V x cj where p and lo are harmonic 
functions. 

Considering the Helmholtz decomposition of 
7 = P7 + Q7, the new term 7 harm can be included 
cither in the divergence-free part or in the curl-free 
part of the decomposition. This decomposition is no 
longer unique, and it is the same for the E/B mode 
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decomposition (24). 

However, the unicity of the reconstruction of ke can 
be obtained under the condition kb = 0. Indeed, if 7 e 
(L 2 ([0, l] 2 ) 2 is split into a divergence-free part P7 and a 
curl-free part Q7, satisfying 



7 = P7 + Q7 and 



ke 




-P7 + Q7, (26) 



then, 7 remains to Te as defined in (12). In addition, the 
convergence ke is uniquely defined, modulo an additive 
constant. Reciprocally, if 7 remains to Te, then we have 
the decomposition 



7 = 




(27) 



and 



K = 



K E 

k b 



-d 2 
di 



d 2 *P + 



(28) 



so kb = — \d\d 2 ip + \d 2 d\ip = 0. The same is true for the 
B modes. Under the constraint ke = 0, 7 remains to Tb 
as defined in (15) and kb is uniquely defined, modulo an 
additive constant. 

This method offers a way to recover the E modes con- 
vergence without the border effects. But the constraint 
kb = should be used carefully, because it assumes that 
there arc no residual systematics in the data due to the ob- 
servational or instrumental effects. We see in the following 
that we can relax this constraint, and impose that k b = 
only on the border, and have a good reconstruction of E 
modes, even in the presence of B modes. 

4. Discretization of the problem 

4.1. Using a harmonic decomposition 

To solve the weak lensing problem on a bounded do- 
main f2 using the formula (26), a possible solution consists 
in writing 7 as 



with 



7 = Vp + Vxw + Vp + Vxu , 



Ap = div7, Aw = curl 7, 



(29) 



(30) 



where p = and co = on the borders of the square SCI 
and 

Ap Q = 0, Auj q = 0, (31) 

i.e., po and luq are harmonic functions on Ct. 

From the relations (24) and (29), we can rewrite the 
E/B modes decomposition as 



K, = 



KE 

k b 



= -(Vp + Vp ) + (V x uj + V x lu ). (32) 



This decomposition is not unique, and we have to find how 
to split the harmonic term 7 harm into Vpo and V x uj . 

By assuming 7 derives from the gravitational potential 
without residual systematics, we can apply to the solution 
the constraint that kb = on all the domain f2. Then, the 
solution to this problem becomes unique and the conver- 
gence is given by Eq. (32). 

Then, po and ujo can be estimated accurately by min- 
imizing the following quantity: 

min ||7- Vp- V x u- Vp - V x uj a \\ 2 + \\k b \\ 2 , (33) 

where Vp and V x oj can be obtained from the relation 
(30) and Vpo an d V x u>q can be obtained by decomposing 
Po and oj in a basis of harmonic functions on (p = 

Z)n=i and w o = J2n=i d n*Pn) and minimizing the 

quantity (33). A set of harmonic functions ^ n can be built 
as 

tp2n-i = n((x 1 +ix 2 ) n ), 

4> 2n = l((x 1 +ix 2 ) n ), for n>l. (34) 

Or, alternatively, 

r/>2n-i = 7e(cxp(n(a;i + ix 2 )/L)), 

tp2n = X(exp(n(a;i + ix 2 )/L)), for neZ*, (35) 

where 1Z and X stand for the real and imaginary parts, 
and L is a constant such that 2-7T L is larger than the size 
of the domain. 



4.2. Using a wavelet decomposition 

To solve the weak lensing problem on a bounded do- 
main f2, another possible solution consists in looking for 
a decomposition of 7 of the form: 

7 = V(p + p )+V x {lu + luq) = Vp' + V xoj'. (36) 
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By considering a zero B modes constraint as in the previ- 
ous section, 

KE = -Vp' + V x uj', with k b = 0, (37) 
J 

the problem can be solved by minimizing the following 
quantity: 

min||7- Vp' - V x lu'\\ 2 + \\n B \\ 2 ■ (38) 

v' : 0J ' 

This minimization problem may involve two bases 
{^jfc } an d {*jfc V } °f the divergence-free and curl-free 
functions on f2 without any boundary conditions: 

V x (c + cjo) = £<Sf V( P + Po ) =EW- 

(39) 

This is the option that has been considered in the pa- 
per. The basis functions 5 , ^" rl and \&4jj. v will consist of 
wavelets. The construction of these wavelets will be pre- 
sented in the next section and in Appendix A. Only a few 
of the boundary elements Wjj? 1 and \&4jj. v do not satisfy 
the conditions Tdivo ■ n = and 7 cur io ' t = on the 
border dil as in the decomposition (30). Then it is possi- 
ble to restrain the zero B mode constraint to these border 
wavelets, and we call this the B-zero border constraint. 

5. Construction of wavelets for Helmholtz 
decomposition 

In this section, we first review briefly the construction 
of divergence-free and curl-free wavelets on an unbounded 
domain. In Appendix A, we present in detail the bases that 
we have selected. Then, we exhibit new divergence-free 
and curl-free wavelets on the square that satisfy suitable 
boundary conditions. 

5.1. Unbounded domain 

Divergence- free and curl-free wavelets 

The compactly supported divergence-free wavelets 
in (L 2 (M. n )) n originally developed by Lemarie-Rieusset 
(1992) are based on the existence of two different ID mul- 
tiresolution analyses (MRA) of L 2 (M.) related by differen- 
tiation and integration. Given a ID multiresolution anal- 
ysis (<pi,tpi) with ipi e C 1+e for e > 0, there is another 
ID MRA (ipo,ipo) such that 

<p' 1 (x) = <p Q (x)-<po(x-l), ip' 1 (x)=4:ip Q (x). (40) 



Different MRA divergence-free wavelets can be con- 
structed from different families of ID MRA's. The con- 
struction of MRA curl-free wavelets follows the same 
approach. Fig. 2 shows an example of a family of ID 
MRA functions and Fig. 3 shows the corresponding 
2D divergence-free wavelets. These wavelets form a ba- 
sis of the space of the divergence-free functions on R 2 . 
Nevertheless, they are not suited to implementing the 
Helmholtz decomposition. 

For this purpose, we prefer to use the hyperbolic 
divergence-free and curl-free wavelets proposed by Deriaz 
& Perrier (2009), which form better conditioned bases. 
They are based on multidimensional hyperbolic wavelets 
introduced by DeVore et al. (1998). The advantage of these 
wavelets is that there is a fast algorithm to compute the 
Helmholtz decomposition in terms of divergence-free and 
curl-free wavelets. 

The Helmholtz decomposition in the wavelet domain 
uses both the divergence-free and curl-free wavelets. As 
the divergence-free and curl-free wavelets form biorthog- 
onal bases in their respective spaces, the coefficients 
dj k rl and dfy corresponding to the divergence-free and 
curl-free parts (eq. (39)) are, in practice, computed with 
an iterative algorithm (see Deriaz & Perrier 2006, for 
more details). 



5.2. Bounded domain 

The Helmholtz decomposition (39) may be generalized 
to bounded domains with non-periodic boundary condi- 
tions by the construction of divergence-free and curl-free 
wavelets on the interval. Since divergence-free and curl- 
free wavelets are constructed from standard, compactly 
supported, biorthogonal wavelet bases, they can incor- 
porate boundary conditions. The following section intro- 
duces a new construction of wavelets well suited to the 
Helmholtz decomposition on bounded domains. 

5.2.1. Wavelet transforms on the interval 

We obtain wavelets on the square by taking tensor 
products of ID wavelets on the interval. The wavelets on 
the interval we introduce here differ from the original con- 
structions (Meyer 1990; Cohen et al. 1993) and rely on the 
lifting scheme (Daubechies & Sweldens 1998). 
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Figure 2. Example of twolD MRA related by differentiation and integration. Scaling functions and associated wavelets for 
biorthogonal MRA with 2 (left) and 3 (right) zero moments. 



tydiv (1,0) 



^div (0,1) 



. t \ 



I TV 



^div (1,1) 



Figure 3. Example of a 2D divergence-free wavelet basis designed from the Lemarie & Rieusset (1992) scheme. 



The properties of these MRA on the interval allow us 
to apply the wavelet Helmholtz algorithm efficiently on 
the square: the partition of the unity property facilitates 
the approximation of the data, and that only ip^ and 
differ from at the boundaries permits us to impose 
boundary conditions on the wavelets. 

An example of spline divergence-free wavelets issued 
from this construction is shown in Fig. 5 for free bound- 
ary (left) and for sliding conditions (right). The bound- 
ary curl-free wavelets can be constructed similarly. We 



applied the resulting boundary divergence-free and curl- 
free wavelets to perform a Helmholtz decomposition on 
the bounded domain. In the next section, we present the 
results obtained with the algorithms that are derived from 
this construction and presented in Appendix B. 



Figure 4. Examples of hyperbolic 2D divergence-free and curl-free wavelets. 
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Figure 5. Example of divergence-free wavelets with a free (left) and a sliding boundary condition (right). 



6. Numerical experiment for weak lensing 

6.1. Simulated data 

We have run realistic simulated convergence mass 
maps derived from N-body cosmological simulations us- 
ing the RAMSES code (Teyssier 2002). The cosmological 
model is chosen to be in concordance with the ACDM 
model. We chose a model with the following parameters 
close to WMAP : tt m = 0.3, <t 8 = 0.9, tt L = 0.7, h = 0.7. 
Each simulation has 256 3 particles with a box size of 162 
Mpc/h. The resulting convergence map covers a field of 2° 
x 2° with 350 x 350 pixels and assumes a galaxy redshift 
of 1. The overdensities correspond to the halos of groups 
and clusters of galaxies. 



6.2. Example with no B modes 

In this experiment, we computed the shear map j s 
from the simulated map n s , and we have extracted a sub 
area j D from 7 s . The lefthand side of Fig. 6 shows the 
simulated k s mass map and its corresponding shear field. 
The righthand side of Fig. 6 shows the subfield r ) D which 
is used for the reconstruction. Then we reconstructed Re 
and kb from -f D , and compared ke to k t , where k t is 
the corresponding original area of n s '. This allows us to 
study how well the reconstructing methods perform with 
the border problem. 

Figure 7 shows the error maps relative to four different 
reconstruction methods: i) the standard FFT approach, ii) 
Seitz & Schneider method, iii) the wavelet Hclmholtz (50 
iterations) with a zero B-mode border constraint, and iv) 
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Figure 6. Left, simulated k s image and its corresponding shear field 7 s . Right, extracted subfield 7°. The extraction breaks 
the periodicity of the data. 




with Fourier Seitz & Schneider 




Figure 7. Reconstruction error maps for four different methods: with the Fourier transform (top left), Seitz & Schneider (top 
right), with the constrained Helmholtz decomposition with the constraint only on the boundary (bottom left), and on all the 
domain (bottom right). 
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the wavelet Hclmholtz (50 iterations) with a zero B-mode 
constraint. The error maps were obtained by e = K rcc — K T . 

We calculated the relative percentage error on both E 
and B modes: eg = ^pr^ ^ an d £b = jj^jf; anc ^ obtain 
e E = 28 %, e £ = 5.3 %, e E = 1.2 % and = 1.0 % for 
the four methods respectively. For the B mode component, 
we obtained = 30.2 % for FFT approach and e B = 1.3 
% for the wavelet Hclmholtz with a zero B modes border 
constraint. Figure 8 shows the reconstructed B-mode map 
for the Fourier method and the wavelet Hclmholtz with a 
zero B modes border constraint. 

This error is not due to the discontinuity on the bound- 
ary because we tested the mirroring technique, which con- 
sists in symmetrizing the data along the (Ox) and (Oy) 
axes in order to recover the continuity in exchange to a 
doubling of the size of the domain. It produced almost the 
same error in terms of norm as a straightforward Fourier 
transform. 

From this experiment, we can conclude that the 
wavelet Hclmholtz method allows us to dramatically re- 
duce the error relative to the border effect, when no B 
modes contaminate our data. 

6.3. Example with B modes 

In this experiment, we have added B modes. The B- 
mode map consists in a Gaussian with a maximum equal 
to 10% of the maximum of k t from the mass map. Since 
the constrained wavelet Helmholtz decomposition assumes 
no B modes, or zero-border B modes, it is interesting to 
see what happens when we introduce some actual B modes 
into the model. 

Figure 9 shows the reconstruction error maps with the 
three methods. We can see that the reconstruction error 
when using Fourier is not amplified in the presence of B 
modes, while the error increases when using the wavelet 
Helmholtz with a constraint on the B modes. This is ex- 
pected since the model is not completely true anymore. We 
note, however, that even if the error increases, it is still 
much below the FFT error. In this case, wavelet Hclmholtz 
with zero border B modes is more robust than wavelet 
Helmholtz with zero B modes. 

Figure 10 shows the original B mode map and the re- 
constructed error maps for both the Fourier transform and 
the wavelet Helmholtz decomposition with a zero B-mode 
border constraint. From this experiment, we can conclude 
that the wavelet Helmholtz method with zero-border B- 



mode constraint is the best trade-off when some residual 
B modes may contaminate the data. 



7. Conclusion 

We have sorted out a direct and original relationship 
between the Helmholtz decomposition and the E/B mode 
reconstruction. The Helmholtz decomposition can be per- 
formed without the FFT, by using divergence-free and 
curl-free wavelets. This wavelet Hclmholtz decomposition 
is iterative, and we have shown that we have different 
strategies to handle the border, thanks to the versatility 
of wavelets. This means that, periodic borders are not a 
curse anymore. Indeed, we were able to design two kinds 
of Helmholtz decompositions, one imposing B modes as 
equal to 0, and another imposing B modes as zero only on 
the border. We have shown that both approaches outper- 
form the FFT-based standard E/B mode reconstruction, 
which suffers from the serious drawback of imposing pe- 
riodic border conditions that concentrate in the low fre- 
quencies. The smaller the field, the more important the 
bias introduced by this periodic border condition. In our 
experiments, using the Helmholtz wavelet decomposition 
reduces dramatically the error. We have shown that the 
wavelet Helmholtz decomposition outperforms also Seitz 
and Schneider method (Seitz & Schneider 1996). Our ap- 
proach also offers the advantage of no needing to pixelate 
the shear map. Indeed, both divergence-free and curl-free 
wavelet transforms can be applied directly on the shear 
catalog without having to pixelate first on a uniform grid. 

However, at this stage, our wavelet Helmholtz cannot 
be properly used yet on real data. Indeed, two main is- 
sues have not been addressed in this paper, the noise and 
the missing data. Both may be problematic: stability, con- 
vergence, and unicity of the solution must be studied in 
detail. This will be done in the future. 
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with Fourier 
"^T 111 - 2 = 30.2% 




Helmholtz (zero-B border constraint) 

"tr"^ 2 = 1.3% 



Figure 8. B-mode reconstructed maps. 




with Fourier 
l|g fT,f = 28.0% 




Helmholtz (B-zero border constraint) 




II 



= 2.1% 



Helmholtz (B-zero constraint) 
'"TT/T"* 2 = 4.0% 



Figure 9. Maps of the reconstruction error on ke for three different methods when B modes contaminate the data. 



Appendix A: Divergence-free wavelets on a square 
domain 

Expression of the divergence-free and curl-free wavelets 

With the ID wavelets of Eq. (40), the 2D isotropic 
divergence-free scaling function is defined as (Lemarie- 
Rieusset 1992) 



$ CUrl (x 1 ,.T 2 ) = 



ipi(x 1 )ip' 1 (x 2 ) 
-ip' 1 (xi)ipi(x 2 ) 



and the details on each scale are obtained from the three 
wavelets 



Horizontal wavelet: 
curl (1,0) f S| ( tpiftxi - ki)(p[(2 j x 2 - k 2 ) 
-V4 (Vxi -fci^i (Vx 2 -k 2 ) 



^^>(x 1 ,x 2 ) 



Vertical wavelet: 
curl (0,1) ln _ ^ _ ( -tpi{2 j xi - ki)ip' 1 (2^x 2 - k 2 ) 



Diagonal wavelet: 



rl(l,l) 



(%1,X 2 ) 



-^' l {Vx l -k l )^ l {Vx 2 -k 2 ) )■ 



The hyperbolic divergence-free vector wavelets are ex- 
pressed as 



(41) 

— Divergence-free wavelets: 



^curl, v _ I 2^(2^X1 - k 1 )M^ J2 X 2 ~ hi) 

j' k ^ X2 > \ -2 jl ip (2 jl xi - ki)ipi(2 j2 x 2 -k 2 ) )' 

(42) 



Curl-free wavelets: 



2^M^ n xi - k!)ipi{2^x 2 - ha) 
2^Vi(2 J1 xi - ki)tp {2^x 2 - k 2 ) 



(43) 
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Figure 10. Original B mode and reconstructed error maps for the B modes map Kg. 



*ij,k(a:i,£2) 
^ij,k(xi,x 2 ) 
*2j,k(a;i,^2) 
*ij,k(a;i,^2) 



(44) 



with the scale parameter j = (ji, j 2 ) G Z 2 and the 
position parameter k = (ki,k 2 ) € Z 2 . The wavelet is 
localized at (2~ : > 1 k 1 ,2~ j2 k 2 ). 

The fast wavelet decomposition in these vector wavelet 
bases proceeds as follows. First, in 2D, we need to define 
two differentiable wavelet bases (Deriaz & Perrier 2006): 

>i(2 J ia;i-fci)V>o(2 j2 a; 2 -fc 2 )' 


M&zi - k x )ip x {2^x 2 - k 2 ) 



M^^i-ki)M^ 2 X2~k 2 ) 


^i(2*si - k 1 )i> {y*x 2 - k 2 ) 

We decompose the field 7 in these bases. Then, as the 

divergence-free and curl-free wavelets are defined by 
^curi = 2,- 3$1 . k _ 2ji 9aJ k) and 

*fk = 2Jl *ij,k + 2-' 2 *2j,k. 

The vector wavelet expansion is obtained by a change of 
basis (cf (Deriaz & Perrier 2006)). 

One can notice that the divergence-free wavelets cor- 
respond to the rotational of the wavelet basis {\{'^ij 1 ,k 1 <& 
ipij 2 ,k 2 )} an d the curl-free wavelets correspond simply to 
its gradient. 

Multi-Resolution Analyses on the interval linked by 
derivation 

We extend this construction to bounded domains using 
special ID wavelets on the interval. To obtain the MRAs 



on [0,1], we modify the filters of multiresolution analyses 
on R on the boundaries. The associated scaling functions 
(father wavelets) form a partition of unity in the sense 
that their sum is equal to one. For instance the quadratic 
spline filter, 

1 3 3 1 

pi(a:) = -<pi(2x+l)+-(p 1 (2x)+-cp 1 (2x-l)+-<p 1 (2x-2), 

(45) 
(46) 



is modified on the left boundary to 



1 



,( 2 ), 



and 



^i 2) (z) = ^i 2) (2x) + j Vl (2x - 1) + -w x (2x - 2). (47) 



^ v~, ■ ^i(2x-l) + -c 

We thus we obtain the corresponding scaling functions 
plotted in Fig. 11. 

To form the divergence-free and curl-free wavelets, we 
need two MRAs linked by derivation. In the following, we 
present a method to integrate an MRA that was partially 
derived from Stevenson (2010). 

We start from a Multi-Resolution Analysis (MRA) 

(<p^\ <Pq \ ■ ■ ■ , ¥>o M ^) 011 t ne interval [0, N], verifying 

M 

"£<p$ ) (x) = l [0iN] (x), and ^ 1) (0)=^ M) (JV) = 1. 



i=l 

We put 



Cq } = j (p { o\x)dx. 
" 



(48) 
(49) 
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Figure 11. Quadratic spline scaling functions adapted to the left boundary of the interval. 



Then we define the MRA (ip^ , <pf \ <p[ M+1) ) by 



^\ x) = l-— I ^>{t)dt 



(50) 



= J (^T)Vt 1] (t) - ^Vo°(*)j (51) 
for i = 2, . . . , M, 



A M+X \x) = ^) fvT\t)dt. (52) 



Then the set {<Pi^} forms a partition of the unity, 

N+n 

^2 tpY\x) = l [0<N] (x), (53) 
i=i 

and satisfies ip^\o) = <p[ M+1) (N) = 1. Thanks to this 
construction, to any wavelet ip{ defined in the MRA 
formed by (ifi , • ■ ■ , <Pi we can associate a wavelet 

ipQ in the MRA formed by (ip^\ ■ ■ . , ^o^) by the relation 

In the example of Section 5.2.1, the corresponding 
MRA (<^o> . . • , ¥>o ) are the linear splines on the inter- 
val. The filter is given by 

<p (x) = ~<po(?x + 1) + M^) + \M2x ~ 1), (54) 
and the left boundary filter by 

$\x)=$\2x) + \<p Qx-l), (55) 



Divergence-free wavelets on a bounded domain 

From the differentiation relation introduced above 

< = 4^, (56) 
we form the boundary divergence-free wavelets: 



*U0,k2)^ X2 ) ~ \ -2^0(2^1)^1 (2 h x 2 - fc 2 



(57) 

For these boundary divergence-free wavelets, the border 
conditions are given by V'i'(O), perpendicular to the border 
and by i/jq(0), tangential to the border. 

Appendix B: Algorithms 

In this appendix, we present a new algorithm based 
on the wavelet Helmholtz decomposition that improves 
the weak lensing E/B mode decomposition. 

Wavelet Helmholtz decomposition algorithm 

The Helmholtz decomposition may be viewed as the 
following orthogonal space splitting: 

(L 2 (R 2 )) 2 = H div0 (K 2 ) ® ± H curl0 (M 2 ), (58) 

where Hdi v o(IR 2 ) is the space of divergence-free vector 
functions and H cur io(I^ 2 ) is the space of curl- free vector 
functions. According to the wavelet decomposition (Deriaz 
& Perrier 2006), we have the following direct sums of MRA 
spaces: 

(£ 2 (M 2 )) 2 =H div0 (M 2 )®H^ iv (M 2 ), 
(L 2 (K 2 )) 2 = H= url (R 2 ) H curl0 (M 2 ). { ™> 

These correspond to the decompositions in the divergence- 
free and curl- free wavelet bases, and, denoting Pdivo as the 
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Figure 12. Linear spline scaling functions adapted to the interval. 



non orthogonal wavelet projector on Hdi v o ~see (Deriaz & 
Perrier 2006) for its definition- we have 



7 = Pdiv 7 + Pdiv 7, 7 = Qcurl 7 



rio7, (60) 



where Pj iv and Ql uri are complementary projectors of 
Pdivo and Q cur i o, respectively. 

Unfortunately, the projectors on the divergence-free 
wavelet and curl-free wavelet bases are biorthogonal pro- 
jectors, so we need an iterative method to provide 7di v o 
and 7 cur i and to write 



7 — 7div0 + 7curl0" 



(61) 



Since the wavelet bases are conditioned well, it is possible 
to apply them recursively in order to construct sequences 
7divo and Tcuric which converge to 7 div0 and 7 curl0 in a 
straight way. Then, in order to decompose the data 7 into 
its divergence-free and curl-free parts 7 = P7 + Q7 , we 

U V 

have implemented the following iterative algorithm: 







and the maximum number of itera- 



1- Set u 1 

tions to I max . 

2- Set i = 0, while H7 — u 1 — v 4 || 2 > e and i < I max then 
iterate: 

1- Set 7ros = 7 - u l - v J 

2- Set u l+1 =u*+P div0 7rcs 

3- Set y rcs = 7 - u' +1 - v 4 

4- Set V l+1 = V'+Q cm ,0 7res 

5- Set i = i + 1 

The algorithm converges when the residual 7* 0S goes 
to and has been proven numerically. We tested this al- 
gorithm with periodic boundary conditions and obtained 
almost the same result as with the Fourier transform. 



Another point is that this algorithm still converges on a 
bounded domain with free boundary condition wavelets 
(see the left panel of fig. ??). 

Weak lensing reconstruction on a bounded domain 

In this section, we present the iterative algorithm that 
aims at adding the constraint kb = in the wavelet 
Hclmholtz decomposition algorithm to improve the recon- 
struction of the E modes convergence in a bounded 
domain: 



7 = P7- 







= -P7 + Q7. (62) 



The algorithm that we have implemented is 

1- Set u° = v° = and the maximum number of itera- 
tions I max . 

2- Set i = 0, while H7 — u 1 — v l || 2 > e and i < I max then 
iterate: 

1- Set 7 rcs = 7 — u l — v* 

2- Set u*+i =u 4 + P div07rc 




3- Set w' = 

4- Set u i+1 = u ! +5 +P di vow 

5- Set 7 ri 4 = 7 - u 4+1 

6- Set v J +5 = v 



v" 

curl Ti'cs 







(-<■)- 









7- Set w H 



8- Set v J+1 = v 1+ 2 - Q cur iow 1+ 3 

9- Set i = i + 1 

This algorithm converges with a slow decay rate. 



16 



Wavelet Helmholtz decomposition for Weak Lensing mass map reconstruction 



An alternative to this algorithm consists in restricting 
the B modes constraint only to the boundary wavelets 
{^fe}- With such an approach, we are able to isolate some 
of the B modes that are not harmonic. The alternative 
algorithm will be 

1- Set u° = v° = and the maximum number of itera- 
tions to I max . 

2- Set i = 0, while H7 — u* — v*|| 2 > e and i < I max then 
iterate: 

1- Set 7 j cs = 7 - u' - v 4 

2- Set U J +2 = U 4 + P d iv0 7rcs with: u * + ^ = 

V 7/ l+5 \l> div 

l~ijk u jk T jk 

3- Compute u£* = £ jfc s ., ^ % ^u^^ 

4- Setw*=( ? + i ) + ( ° ) 

V - u su 2 J V v sn 2 J 

5- Set u l+1 = u i+ 5 + Pdivo w J 

6- Set 7^/ = 7 - u 4+1 - v J 

7- Set v l +3 = v* + Q curio7 ri' with v i+ 5 = 

8- Compute v^ 5 = £\ fc s t vj**^ 1 

9- Set w*+^ = C ° ) + ( ilk) 

10- Set V l+1 = V l +5 - Qcurl W i+ 5 
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